#生成测试数据
#4个状态，2种结果
import numpy
from hmmlearn import hmm
#模型训练
# A=numpy.array([[0,1,0,0],[0.4,0,0.6,0],[0,0.4,0,0.6],[0,0,0.5,0.5]])
# B=numpy.array([[0.5,0.5],[0.3,0.7],[0.6,0.4],[0.8,0.2]])
# pi=numpy.array([0.4,0.2,0.3,0.1])
# t=5#规定t=10
# n=10000#规定数据量
# s=numpy.empty((n,t),dtype=int)
# v=numpy.empty((n,t),dtype=int)
# for i in range(n):
#     #生成初始状态
#     start_s=numpy.sum(pi.cumsum()<numpy.random.rand(1))
#     #根据初始状态生成其他状态和对应观测值
#     s[i,0]=start_s
#     v[i, 0] = numpy.sum(B[start_s].cumsum()<numpy.random.rand(1))
#     for j in range(1,t):
#         start_s = numpy.sum(A[start_s].cumsum() < numpy.random.rand(1))
#         s[i,j]=start_s
#         v[i, j] = numpy.sum(B[start_s].cumsum()<numpy.random.rand(1))
# model=hmm.CategoricalHMM(n_components=4, tol=1e-10,n_features=2,startprob_prior=numpy.ones(4)/4,emissionprob_prior=numpy.ones((4,2))*0.5,transmat_prior=numpy.ones((4,4))/4)
# model.fit(v.reshape(-1,1),lengths=numpy.ones(n,dtype=int)*t)
# print(model.transmat_)
# print(model.emissionprob_)
# print(model.startprob_)
# predict_s=model.predict(v.reshape(-1,1),lengths=numpy.ones(n,dtype=int)*t)
# print(predict_s)
A=numpy.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])
B=numpy.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])
pi=numpy.array([0.2,0.4,0.4])
model=hmm.CategoricalHMM(n_components=3,n_features=2)
model.transmat_=A
model.emissionprob_=B
model.startprob_=pi
X=numpy.array([[0],[1],[0]])
#概率计算
print(numpy.exp(model.score(X,lengths=numpy.array([3]))))
#解码
print(model.predict(X,lengths=numpy.array([3])))

#公式推导参考：
# https://blog.csdn.net/Elenstone/article/details/104902120
# https://zhuanlan.zhihu.com/p/85454896       参数估计时要考虑多个序列样本，不仅仅是一个
